function [MeshData] = ReadGmshTA(file_name)

fid = fopen(file_name);

tline = fgetl(fid);

while(~strcmp(tline,'$Nodes'))
    tline = fgetl(fid);    
end %while

tline = fgetl(fid);
numNodes = str2double(tline);
tline = fgetl(fid);

x = zeros(numNodes,1);
y = zeros(numNodes,1);
z = zeros(numNodes,1);
nodeIDs = zeros(numNodes,1);
index = 1;

while (~strcmp(tline,'$EndNodes'))
    tmpv = sscanf(tline,'%d %f %f %f');
    nodeIDs(index) = tmpv(1);
    x(index) = tmpv(2);
    y(index) = tmpv(3);
    z(index) = tmpv(4);    
    index = index + 1;
    tline = fgetl(fid);
end %while

while(~strcmp(tline,'$Elements'))
    tline = fgetl(fid);    
end %while

tline = fgetl(fid);
totNumEle = str2double(tline);
tline = fgetl(fid);

eleNodes = zeros(totNumEle, 4);
eleNodes2D = zeros(totNumEle, 3);
bdNodes = zeros(totNumEle,1);
phys = zeros(totNumEle,1);
phys2D = zeros(totNumEle,1);

num3DEle = 0;
num2DEle = 0;
numBDnodes = 1;  %remember that BDnode count is 1 greater than actual found

while(~strcmp(tline,'$EndElements'))
    tmpv = sscanf(tline,'%d');
        
    if(tmpv(2) == 2)  %if 2D element
        num2DEle = num2DEle + 1;        
        for ii = 2:-1:0
            if(~sum(bdNodes(1:numBDnodes) == tmpv(end-ii)))                
                bdNodes(numBDnodes) = tmpv(end-ii);
                numBDnodes = numBDnodes + 1;
            end %end if
        end %for ii
        
        eleNodes2D(num2DEle,:) = tmpv([end-2, end-1, end]);
        phys2D(num2DEle) = tmpv(4);
        
    elseif(tmpv(2) == 4)  %if 3D element
        num3DEle = num3DEle + 1;
        eleNodes(num3DEle,:) = tmpv([end-3, end-2, end-1, end]);
        phys(num3DEle) = tmpv(4);
    end %if   
    
    tline = fgetl(fid);
end %while
fclose(fid);

if((num3DEle + num2DEle) == totNumEle) 
    display('Sanity Check: Number of Elements... ok');
else
    display('Sanity Check Failed: Number of 2D/3D elements is inconsistent');
end %end if

% Saving Mesh Information to output structure
MeshData.nNodes = numNodes;
% MeshData.BdNodes  = bdNodes(1:(numBDnodes-1), 1);
MeshData.xNodes = x;
MeshData.yNodes = y;
MeshData.zNodes = z;
% MeshData.NodesID = nodeIDs;
MeshData.nElements3D = num3DEle;
MeshData.nElements2D = num2DEle;
MeshData.EleMatrix3D = eleNodes(1:num3DEle, :);
MeshData.EleMatrix2D = eleNodes2D(1:num2DEle, :);
MeshData.PhysNumber3D = phys(1:num3DEle, 1);
MeshData.PhysNumber2D = phys2D(1:num2DEle, 1);

%also store centroid of each element in (x,y,z) format
MeshData.centroids = zeros(MeshData.nElements2D, 2);
for n = 1:MeshData.nElements2D
    MeshData.centroids(n, 1) = sum(MeshData.xNodes(MeshData.EleMatrix2D(n,:))) / 3;
    MeshData.centroids(n, 2) = sum(MeshData.yNodes(MeshData.EleMatrix2D(n,:))) / 3;
end %ii   


